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ABSTRACT 

Observations of dusty debris disks can be used to test theories of planetesimal coagulation. Plan- 
etesimals of sizes up to a couple thousand kms are embedded in these disks and their mutual collisions 
generate the small dust grains that are observed. The dust luminosities, when combined with in- 
formation on the dust spatial extent and the system age, can be used to infer initial masses in the 
planetesimal belts. Carrying out such a procedure for a sample of debris disks around Sun-like stars, 
we reach the following two conclusions. First, if we assume that colliding planetesimals satisfy a 
primordial size spectrum of the form dn/ds oc s~'^, observed disks strongly favor a value of q between 
3.5 and 4, while both current theoretical expectations and statistics of Kuiper belt objects favor a 
somewhat larger value. Second, number densities of planetesimals are two to three orders of magni- 
tude higher in detected disks than in the Kuiper belt, for comparably-sized objects. This is a surprise 
for the coagulation models. It would require a similar increase in the solid surface density of the 
primordial disk over that of the Minimum Mass Solar Nebula, which is unreasonable. Both of our 
conclusions are driven by the need to explain the presence of bright debris disks at a few Gyrs of age. 
Subject headings: circumstellar matter — Kuiper belt — planetary systems: formation 



1. INTRODUCTION 

Dusty disks made up of rocky and icy debris have been 
obser ved around other sta rs, both in reflected optical 
light (jSmith fc Terrilelll984f) and in lo ng wavelength ther- 
mal radiation (|Aumann et al .1119841 ). Multiple surveys 
have reported that a significant fraction of main-sequence 
stars harbo r detectable infrared excesses: ~ 15% solar- 
type stars (iTrilling et al .' 2008; ' Lawler et all 120091 ) . and 
- 30% for A-stars (jSueLaL 2QQ1) . The infrared lumi- 
nosity, when compared to the luminosity of the central 
star, ranges from ~ 10~^ to ~ 10"'^. In contrast, the 
fractional dust lumino sity from the Kuipe r belt is esti- 
mated to be - 10""^ (|Teplitz et al.lll999D and remams 
undetected. 

The observed excess luminosities arise primarily from 
small (~ p, m — mm) dust grains. Due to t heir short sur- 
vival time ()Artvmowicz &: ClampinI [l997[ ). these grains 
are believed to be continuously produced by collisions 
between large parent bodies ('planetesimals'). These 
planetesimals, analogous to the Kuiper belt objects in 
our own system, are in turn left-overs from the epoch of 
planet formation. 

In this article, we describe how we can use debris disks 
to test theories of planetesimal formation. We first fo- 
cus our attention on the primordial size spectrum of 
planetesimals, often characterized by a single power-law, 
dn/ds oc s~'^, where s is the size. In the following, we 
briefly summarize theoretical understandings and obser- 
vational evidences for the value of q. 

The conventional picture of planetesimal formation is 
composed of a number of steps. The formation of the 
first generation planetesimals is not yet well- understood 
and is an area of active res earch (see, e. g.lYoudin fc Shu 
200aiDominik et al.ll20'07t iJohansen et al.ll2007HGaraud 
2007[) . If these are sufficient ly massive, gravity dominate s 
their subsequent growth (jWeidenschilling et al.l I1997D . 
At first, objects grow in an orderly fashion, where col- 



lisions and conglomerations occur at rates that are pro- 
portional to their geometric cross sections. But when 
these bodies become so massive that the effect of grav- 
itational focusing becomes significant, run-away growth 
commences where the largest bodies accrete small plan- 
etesimals at the highest rate and quickly distance them- 
selves from their form er peers (iWetherill fc Stewarti[l98l 
iKokubo fc Idal [19961) . The run-away phase is succeeded 
by the oligarchic phase where individual large bodies are 
responsible for stirring the small bodies that they accrete 
(jKokubo fc Ida 1995. 1998). At the end of these steps, 
an entire size spectrum of planetesimals are produced. 
This is the 'primordial spectrum'. 

During the run-away phase, N-b ody simulations have 
typically produced a slope of g ^ 6 (jKokubo fc Idal[l996l : 
[Morishima et al. 20]^. This slope is naturally explained 
if there is ener gy equi-partition am ong planetesimals of 
different sizes (jMakino et al.l[T998| ). Moreover, one ex- 
pects that the distribution becomes shallower (smaller q) 
if larger planetesimals have higher kinetic energies. This 
indeed occurs during the oligarchic phase when all small 
and intermediate-sized planetesimals are stirred to the 
same v elocity dispersion. The value of q is then reduced 
to w 4 (|Morishima et al.ll2008| ). 

Using par ticles-in-a-box s i mulat io ns and later hybrid 
simulat ions. iKenvon fc Luul (|1999f ): iKenvon fc BromlevI 
()2004bll2008l) followed the growth of planetesimals. They 
also found that q decreases with time after the run- 
away phase, finishing up with 3.75 < q < 4.5 for 
planet esimals of sizes be t ween 10 and 1000 kms. Ke- 
centlv. lSchlichting fc Saril()2011l ) argued analytically that 
a q = 4 spectrum is the natural outcome of conglomera- 
tion. 

Observational constraints on the value of q currently 
come exclusively from counting large Kuiper belt ob- 
jects. Kuiper belt objects larger than about 30 — 50 kms 
are commonly believed to be primordial. Collision 
timescales for these bodies well exceed that of the So- 
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lar s ystem age (jPavis fc Farinellal 119971 : iBianco et alj 
IMoh . The size distribution for these bodies can 
be probed by present-day surveys. Pubhshed val- 



ues for q are scattered: q = 4.0_ 
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4.25 ± 0.25 (|Fraser et all l20M) . 
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4.5 ± 0.4 (iFrascr fc Kavelaad Hool) and q = 4.5: 
(|Fuentes fc Hohnan 200^ . This scatter may be intrinsic 
and reflect both the different size ranges and the differ- 
ent d ynamical populations emphasized by various sur- 
veys (jSernstein et al.ll2004l : lDonnisonll2006l : iFraser et al.l 
|2010( ). For bodies smaller than ~ 30 kms, the size dis- 
tribu t ion adopts a shallower power-law ([B ernstein ct afj 
[200l iFuentes fc HolmanI [2008t ISchhchting et al.l 120091) . 
This break in the power-l aw index has been argued to be 
due to coUisional erosion (IPan fc Sarill2005D. but a diff er- 
ent opinion has surfaced ( Charnoz fc Morbidellil[2C)06( ). 

So at least for the value of current coagulation mod- 
els appear to be vindicated by the observations. These 
models enjoy a further success. In the Kuiper belt re- 
gion, the solid mass of t he so-called Minimum Mass So- 
lar N ebula is ~ 10 (jHavashil [T98lt IWeidenschillind 
Il977| ) , while the mass in large Kuiper belt objects is 
estiraated to be < 0.1 M m (see, e.g. iGladman et al.l 
120011: iBernstein et al.l[200l . This large difference, how- 
ever, is explained by current models where the forma- 
tion of large planetesimals has a very low effic iency 
(|Bromlev fc Kenvonl[200llSchhchting fc SarillMl . 

With these two remarkable concordances, one wonders 
if debris disks will ever tell us anything new and unex- 
pected. Furthermore, every debris disk likely has a differ- 
ent initial condition and evolves in a different dynamical 
environment. For instance, dynamical interactions with 
Neptune or other planets may ha ve qualitatively affec ted 
the evolution of the Kuiper belt (jLevison et al.ll2068l ) . It 
seems difficult, therefore, to extract any universal truth 
about the formation process from these disparate objects. 

However, based only on a modest sample of debris 
disks, we argue in this paper that there is already a seri- 
ous issue in current coagulation models. 

To achieve this, we first construct a simple coUi- 
sional model ([J2]) to compare against the set of de- 
bris disks reported in iHillenbrand et aH ()2008f ). Our 
collisional mo del does not differ i n essence from pre- 
vious works (iKrivov et al.l 120051: IWvatt et all l2007bl: 
iLohne et al.ll2008() . but we interpret the observations in a 
new way. This allows us to measure the value of q as well 
as the initial masses of planetesimal belts (|j4]). The lat- 
ter result challenges the current models of planetesimal 
formation ([S]). We summarize in ^ 

2. MODEL: LUMINOSITY EVOLUTION OF A DEBRIS DISK 

The debris phase commences when eccentricities of the 
primordial planetesimals are further increased so that 
they no longer coalesce at encounter, but are instead bro- 
ken into fragmentsQ In this phase, the smallest primor- 
dial planetesimals enter into a collisional cascade first, 
followed by progressively larger bodies. During the col- 
lisional cascade, a primordial body is broken down into 
smaller and smaller fragments until all its mass ends up 
in small grains. The small grains may spiral in towards 
the star due to Poynting-Robertson drag, as happens in 



^ IKenvon fc BromlevI l|2008l ) find that fragmentation begins once 
Pluto-sized bodies form. 



the Solar system, or, be ground down by frequent colli- 
sions to sizes so small that they are promptly removed 
by radiation pressure, as happens in bright debris disks 
(|WvattJl2005h . 

2.1. Debris Rings 

We model the debris disk as a single, azimuthally 
smooth ring composed of planetesimals of different sizes. 
The ring is centered at a semi-major axis a with a full 
radial width of Aa and a constant surface density. We 
take Aa/a = 0.1 ^ 1 as our standard input. This is mo- 
tivated by the following observations. Spatially resolved 
debris disks often appear as narrow rings. Examples are, 
Aa/a ^O.l for AU Microscopii dF itzgerald et al. 200^^, 
- 0.5 f or HP 10647 (ILiseau et al.l (20100. - 0.3 for HD 
92945 (iGolimowski et al.l l2007D. ~ 0.3 for HD 139664 
(IKalas et al.ll2006D . ^ 0.2 for HD 207129. (IKrist et al.l 
120101) , ^ 0. 5 for e Eridani (IDent et aLl[2000l) . y 0.1 for 
Foma lhaut ()Kalas et al.l 120051 ). ^ O."2lor Vega (jSu et all 
120051 ). Similarly, unresolved disks often exhibit spectral 
energy distr ibution that is well fit by a single temperature 
blackbodv jHiUenbrand et al.l l2008t iNilsson et all I2OIOI: 
iMoor et al.ll2011l) . This ring- like topology also show up 
in our own Solar system, hence the name the asteroid 
"belt" and the Kuiper "beh" . 

2.2. Initial Size Distribution of the Planetesimals 

We adopt the following power-law forms for the initial 
size distributions. 



dn 

ds 



-«3 



(1) 



The index q^ is the primordial size index for large bod- 
ies, like one that arises out of conglomeration models. 
Previous studies of collisional debris disks have taken 
this value to be a given, in fact it is commonly set 
to be th e power law one expects f rom collisiona,l equi- 
librium (IKrivov et al.l l2005l I2OO6I: IWvatt et al.l I2007bl : 
ILohne et al.l 120081 ). In contrast, in this contribution we 
use the observed sample to measure this value. 

In equation ([T]), Sbig is the size of the biggest plan- 
etesimals, Sinin the smallest. The intermediate size Sgmaii 
is introduced for the purpose of mass accounting: the 
original mass counts only those between Sbig and Ssmaii, 
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While Smin is naturally taken to be the size at which 
radiation pressure unbinds dust grains from the star ( 
^ /.jm for a Sun- like star), we discuss our choice for Sbig 
and Ssmaii below. 

Motivated by the observational and numerical results 
discussed in ^ we investigate values of q^ between 3.5 
and 5. The value 53 = 4 has the special property that 
mass is distributed equally among all logarithmic size 
ranges, while masses in systems with q^ > 4 diverge to- 
ward the small end. The intermediate size Ssmaii is intro- 
duced, partly to avoid dealing with this divergence. For 
sizes below Ssmaii, we assume that collisions have set up 
an equilibrium power law with index qi (see Appendix) . 
So, the intermediate size Ssmaii can also be interpreted as 
the collisional break size at time zero. For our study, we 
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set Ssniaii — 100 m. For our typical disks, wc find that, 
within a few million years, collisional equilibrium is es- 
tablished for bodies up to sizes ~ 1 km. So the choice of 
Ssmaii is not important for late time evolution. 

The choice of size for the largest bodies, Sbig, deserves 
some discussion, as it affects the qualitative character of 
the evolution. As a collisional cascade progresses, bod- 
ies of larger and larger sizes come into collisional equilib- 
rium, opening up fresh mass reserve to produce the small 
particles. Once the largest bodies enter into collisional 
equilibrium, t he dust production rate decays with time 
as Ltr oc IW vatt c t al.l[2'007al). Two previous stud- 
ies (jWvatt et al..,2007bnL61me et al.ll2008D have adopted 
sizes for the largest bodies of Sbig = 30 and 74 km, re- 
spectively. For some of their disks, the largest bodies 
can enter collision equilibrium during the lifetime of the 
system. 

Both Kuiper belt observations and numerical studies of 
coagulation favor a largest size of ^ 1000 km. The largest 
object yet found in the Kui per Belt, (136199) Eris, has 
a radius of 12 00 ± 50 km ferown et al. 2009). In the 
simulations of IKenvon fc BromlevI (I2004ai) . coagulation 
of planetcsimals at 30 - 150 AU produces bodies as large 
as 1000 - 3000 km. When the largest bodies reach this 
size, self-stirring increases the velocity dispersion and col- 
lisions become destructive rather than conglomerating. 

Therefore, we adopt a maximum body size of 1000 km 
in our study. Our quoted masses reflect this choice of 
Sbig- Our largest bodies never enter into collisional equi- 
librium. If this assumption turns out to be erroneous, 
namely, Sbig is much smaller and enters into collisional 
cascade within system lifetime, our model would under- 
estimate the initial masses for old disks. As a result, we 
would overestimate the value for (73. 

2.3. Collisions 

We only consider collisions that are catastrophically 
destructive. A catastrophic collision is defined as one 
that removes at least 50% of the mass of the primary 
body. In so doing, we have implicitly assumed that 
both cratering collisions and conglomerating collisions 
are unimportant. When a destructive collision occurs, 
the total mass (bullet plus target) is redistributed to all 
smaller sizes according to dn/ds oc s~^. This choice is 
somewhat arbitrary and we have confirmed that modi- 
fying it (within reasonable bounds) does not change our 
results. 

We do not model evolution of the orbital dynamics as 
bodies collide. This is justified by the discussions in H5.2[ 

Let the chance of collisions between two bodies of sizes 
s and s' be, 

J collision — „ a , i 1*^7 

Here, 27raAa is the surface area spanned by the debris 
ring in the orbital plane, and torb is the orbital period. 
Gravitational focusing is negligible for the high random 
velocities we consider here. The typical encounter veloc- 
ity, for particles with eccen tricity e and inclination i, is 
(|Wetherill fc StewartJUOOl 

Wcoi = \/l.25e2-hz2wkop, (4) 
where Vkep is the local Keplerian velocity. We adopt 



i ~ e/2 so Dcoi ~ l-32eukGp- As argued in i i5.21 it is 
reasonable to assume a constant eccentricity (and incli- 
nation) for all bodies. We take a value of e ~ 0.1 as the 
standard input, and discuss this assumption in fj5l 

We denote the specific impact energy required to catas- 
trophically disrupt a body (target) as Q*. The scaling 
of Q* with the size of the target depends on whether 
its strength is dominated by material cohesion or self- 
gravi ty. We adopt the following form ()Benz fc Asphaud 
|l99l), 

Q* = A(-^Y + Bp(-^Y (5) 
V 1 cm / VI cm / 

where p is the bulk density which we take to be 2.5g/ cm"^. 
The first term on the right-hand-side describes the inter- 
nal strength limit, important for small bodies, while the 
second term the self-gravity limit, important for larger 
bodies. 

The strength law sets the size of the smallest bullets 
required to destroy a target. Since these are also the most 
numerous, they determine the downward conversion rate 
of mass during a collisional cascade. As such, the power 
indexes in the strength law directly determine the size 
spectrum at collisional equilibrium. For a strength law 
of the form Q* oc s'^, the equilibrium size spectrum is 
dn/ds oc s"9, with (jPurda fc Dermottl[T997[ ): 

g=(21 + c)/(6 + c). (6) 

The famous Dohnanyi-law (jPohnanvil Il969[) . dn/ds oc 
s~^'^, obtains from c = 0. 

The value and form for Q* are notoriously difficult to 
assess. It depends on. among other factors, material 
composition, porosity and impact velocity. A number 
of computations and compilations have appeared in the 
literature. We select three representative formulations 
for our study (Fig. [ij. 

Based on a variety o f exper imental data a nd SP H sim- 
ulations, iKrivov et all (|2005f ): lL5hn e et aD (|2008[) advo- 
cated the following choices, A = 2x lO^erg/g, a = —0.3, 
B = 0.158, /3 = 1.5. We call this the 'hard' strength law. 
In this case, the collision spectrum satisfies q « 3.6 and 
3.0, in the strength and gravity r egimes re spectively. 

Based on energy conservation, iPan fc Sar i (,20051 ) cal- 
culated a destruction threshold for bodies that have zero 
internal strength and obtained B = 3.3 x 10~^, /3 = 2. 
So bodies at 100 km is we aker by a fac t or ~ 1000 than 
their counterparts in the IKrivov et al.l (j2005D formula- 
tion. Wc refer to this as the 'soft' strength law. A softer 
strength implies smaller bullets and therefor e more fre- 
quent destruction of the targets. iPan fc Saril (j2005f ) did 
not consider smaller bodies that are strength bound. We 
adopt A = 2 X 10''erg/g and a = —0.3 in this range to 
c omplete the soft presc r iption . 

iStewart fc Leinhardtl (|2009f ) proposed a strength law 
that depends on impact velocity, 

Q* = (500s-"-33 + 10-^1-2) i;0„f, (7) 

For a typical velocity Vco\ = 500m/s and for bodies 
greater than 1km, this gives rise to a strength law that 
falls in-between that of the hard and the soft case. 
We call this the medium strength law. Note that this 
strength law is much weaker than the other two for small 
bodies. 
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radius (cm) 

Fig . 1 — Prescr i ption s for specific strength from ILohne et al.l 
|[200g'l . IPan fc Saril II2005I') and IStewart fc LeinhaTatl (WOS) . plot- 
ted here as functions of target sizes. We insert an impact velocity 
of 500 m/ s to evaluate the last prescription. Strength of small bod- 
ies are dominated by material cohesion, while that of larger bodies 
by self-gravity. Transitions between the two limits occur around 
100 m (the hard and the medium laws) or around 10 km (the soft 
law) . Strength for bodies smaller than 1 cm are extrapolations as 
both laboratory and numerical experiments only concern bodies of 
larger sizes. 




system age (yrs) 

Fig. 2. — Time evolution of the break-size in a model system, 
with Mo = 12. 4M©, qs = 4.0, e = 0.1, a = 31AU, and Aa/a = 0.1. 
Here, break-size is defined as the size at which all bodies initially at 
that size have encountered of order one destructive collision. Break 
size increases with time monotonically as larger bodies enter into 
coUisional cascade. The numerical results are shown as solid curves, 
while the analytical scaling relations (see Appendix) are plotted as 
dashed lines. The bends in the curves occur at s ^ si, i.e., sizes 
for which material cohesion and gravity binding are comparable. 
The set of thick curves arc for the case of hard material strength, 
while the thin lines for soft strength. 

For the strength laws we consider, transitions from ma- 
terial strength domination to self-gravity domination oc- 
cur at size s « si, with si ranging between 100 m (the 
hard and the medium laws) and 10 km (the soft law). 

2.4. Luminosity Evolution 

The planetesimal disk, starting from an initial disk 
mass of Mq, and an initial size spectrum (eq. [IJ, is nu- 
merically collided and ground down. We divide the par- 
ticles between Sgmaii and Sbig into 500 equal logarithmic 
size bins. The time-step for the simulations is adaptively 




system age (yrs) 

Fig. 3. — Evolution of fractional luminosity, Lir/L*, for the 
system in Fig. [2] The thick line is obtained using the hard strength 
law, and the thin line the soft one. The evolution proceeds in two 
stages: the flatter early stage when coUisional cascade only involves 
small bodies that are bound by material cohesion; and a steeper 
later stage where bodies bound by self-gravity enter the cascade. 
At late times, fractional luminosity decays as t""-^ (eq. ||A9||) . 
While break-sizes differ for the two adopted strength laws (Fig. [2| , 
this appears to have little influence on the overall luminosity. 

set SO that over one time-step, the maximum mass gain 
(from larger bodies) or loss (to smaller bodies) per bin 
falls below 5%. The net mass change is substantially 
smaller than this due to the cancellation between gain 
and loss. 

We calculate the fractional brightness of the dust disk, 
Lir/ , by integrating the geometrical cross section over 
all grains. This assumes that grains are perfect absorbers 
at the optical and can emit efficiently in the infrared. 

An example of such a calculation is reported in Figs. 
[2] & [31 To understand these results, a simple analytical 
model (see Appendix) is introduced. Scaling relations 
obtained using this analytical model compares well with 
our numerical results. 

Fig. [2] shows that, with time, larger and larger plan- 
etesimals enter into coUisional cascade. Within a million 
years or so, the cascade has advanced to size of order 
one kilometer. Beyond this time, bodies bound by self- 
gravity can be gradually eroded. By 1 Gyrs, bodies with 
sizes 10 — 100 kms may be affected. The exact value 
depends on the strength law. The dust luminosity is re- 
lated to the dust mass, which is in turn related to the 
dust production rate. The dust production rate, on the 
other hand, is simply the primordial mass stored at the 
break-size divided by the system age. If the primordial 
spectrum is such that a large amount of mass is piled at 
the large end, debris disks would not exhibit significant 
fading even up to a few billion years. 

Fig. [3] shows that dust luminosit y Lm /L^, oc 
for 53 = 4, consistent with equation (|A11[) . That same 
equation also demonstrates that the value of B, strength 
constant for bodies bound by self-gravity, affects the lu- 
minosity only minorly. This is born out by results shown 
in Fig. O 

An important result on which we base our later anal- 
ysis is shown in Fig. Luminosity evolution for disks 
with the same initial mass but different are depicted. 
As eq. ((X9|) predicts, L cx t(93-3)/(2-g3)_ jf jg shallow 
(e.g. qa < 4), most of the initial mass is deposited at the 



5 



10"' 



10' 



q3=3.5 
q3=4.0 
q3=4.5 
q3=5.0 



1 00 Myrs 

system age (yrs) 



1 Gyrs 



Fig. 4. — Luminosity evolution for disks witli different gs but 
the same initial mass (SMq). Systems with a steeper primordial 
size spectrum (larger 53) exhibit a more pronounced decline of 
luminosity with time, since at a given time, a bigger fraction of 
their mass reservoir has been depleted. Systems with shallower 
gs (e.g., ga = 3.5), on the other hand, are initially dimmer due 
to the relative shortage of smaller rocks, but eventually outshine 
the higher 93 disks as they can hold on to their mass reservoir for 
longer. The luminosity decay of observed disks that span a large 
range of ages can thus be used to infer the value of 53. All other 
parameters here are similar to those used in Fig. [3] and we adopt 
the hard strength law. 

largest planetesimals. This mass reservoir is harder to 
reach by coUision and allows the disk to remain brighter 
at later times. In comparison, disks with a steeper 53 
decay faster. 

If one observes a collection of debris disk all at the 
same age, intrinsic scatter in, e.g., initial masses, makes it 
impossible to differentiate between models of different q^. 
However, a collection of disks with a large age spread can 
be used to constrain ^3. This we proceed to demonstrate. 

3. OBSERVED ENSEMBLE 

Severa l debris disks surveys have been carried out 
see, e.g. iSu et al.ll2006HTrilling et al.|[200lirawler et al.l 



120091: iMoor et all 120 11^" . The sa mple of most in- 
terest to us is that reported in iHillenbrand et al 



( 2008 ). Togethe r with upda tes in iCarpenter et al 



2009( ). IffiUcnbran d et al.l (|2G08[ ) presented a collection 
of debris disks around F/G/K type stars, obtained as 
part of the Spitzer program on Formation and Evolution 
of Planetary Systems (FEPS). This sample is unique in 
that both the stellar age and the radial distance of the 
dust ring are determined: isochrone fitting provides the 
age for the host stars (spanning from ~ 10^ years to a 
few 10^ years), while multi-band photometry and spec- 
tral energy fitting yield the semi-major axis of the dust 
ring. Together with fractional luminosity of the dust 
belt, these provide the most important constraints to in- 
fer the primordial properties of parent planetesimals. 

To obtain the blow-out size (smin) for each system, we 
take luminosity values f or the central stars as given in 
(jHillenbrand et all 120081 ) , and we assign stellar masses 

1 /3 

by assuming that M^, cx , as appropriate for solar 
type main-sequence stars. 

Out of the 31 disks listed in IHillenbrand et al.l ()2008D . 
we focus only on a sub-sample of 13 disks that appear 
radially uncxtcndcd and arc around main-sequence stars. 
In .Hillenbrand et al.. (,2008i ) , emission from each disk is 
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Fig. 5. — Cumulative distribution of stellar ag e (left pane l) and 
dust luminosity (right panel) for the Hillenbrand et al.l l|2008l ^ sam- 
ple. The extended systems (solid curve) tend to be younger and 
brighter than the unextended systems (dashed curves). 

initially fitted with a single temperature blackbody (a 
ring). If agreement between the 24/im/33/im fit and the 
3Sfim/70fim fit is poor, they argue that the disk is likely 
radially extended and fit the data instead with two radial 
components. Since our numerical model is a one-zone 
model, we find that including the extended sources into 
our analysis causes significant scatter in the results. This 
leads us to discard them for the current analysis. We 
have excluded HD 191089 from our sample. Its fluxes 
in 13 /im and 33 fim arc not measured, and cannot be 
reliably identified as an unextended source. In all, we are 
left with 13 sources. 

It is interesting to note that most of the extended 
sources are relatively young, all younger than a few hun- 
dred million years. In contrast, the unextended sources 
have a larger age spread, lasting till a few billion years 
(Fig. [5]). All systems may be born with more than one 
debris rings, but after a sufficiently long time, only the 
outermost ring, which has the longest erosion timescale, 
remains shining. The extended system are also brighter 
than the average, likely related to their relative youth. 

4. THE PRIMORDIAL SIZE SPECTRUM REVEALED 

We have a simple strategy. Knowing the luminosity, 
the age and the semi-major axis of each debris ring, we 
use our coUisional model to backtrack the evolution to 
infer its initial mass in the planetesimal belt. These ini- 
tial masses, when plotted against system ages, should 
show a spread. One expects this spread to be constant 
across all ages, as disks formed at different cosmic times 
likely have the same distribution of disk properties. This 
property could be used to test model assumptions. How- 
ever, using the spread is difficult due to selection effects. 
For instance, low mass disks may become too dim at late 
times to be observable. So we propose instead to study 
the upper envelope of this spread. The upper envelope 
should be flat with age for the correct model. From our 
analytical scaling relations (see Appendix), we find that 
the most important parameter in our model that affects 
this mass slope is ga, the power- law index in the primor- 
dial size spectrum. 

The results of such a procedure are shown in Fig. [HI 
Models with = 4.5 or greater appear to be excluded 
by data, as they would require a rise of initial disk mass 
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Fig. 6. — Inferred disk initial masses, plotted against system ages, for the uncxtended systems in [Hillenbrand et al.l 1)20081 ). The four 
panels present four different choices of 93. The other parameters chosen are e = 0.1, Ssmall = 10*cm, and Aa/a = 0.1. A value of 
93 G [3.5,4.0] is preferred: the upper envelopes for the disk mass remain constant at all ages in the two top plots. Models with higher 53 
are excluded as they require a rising upper envelope. In addition, the 93 = 5 model requires unphysically large disk masses for very old 
disks. 
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Fig. 7. — Similar to Fig. [B] except 53 is fixed at 4 and a number of parameters are varied to test how the inferred initial masses depend 
on them. Specific values for the inferred mass may change, when the small end of the primordial size spectrum (Samalli top-left), the 
eccentricity of particles (bottom- left), the fractional width of the debris ring (top-right), and the adopted strength law (bottom-right) are 
varied for all systems. However, the important indicator for our study, the upper envelope of the masses as a function of system age, 
remains fiat. So the conclusion that 53 ~ 4 remains valid. 



7 



10" 



^ 10-= 

tn 
(0 
E 

« 10 



= 1-10km:q3=7; 10-10^ km: q3=4 

- A 1-10km: q3=5.5; 10-10^ km: q3=4 O 



o 



o 



10" 



100 Myrs 1 Gyrs 

system age (yrs) 

Fig. 8. — Inferred initial masses for a broken power-law size- 
distribution. We investigate t wo particular forms , mot ivated 
by coagulation simulati ons by IKenvon fc BromlevI I I2008I ) and 
ISchlichting fc Saril 11201 II ). respectively. Other parameters adopted 
are e = 0.1, Aa/a = 0.1, and the hard strength law. The inferred 
disk mass rises sharply with system age. Moreover, to make old 
and bright systems, we require disk masses that approach the mass 
of Jupiter. 

with stellar ages. The reason behind this is transparent 
by studying Fig. |4l Models with ~ 3.5 and 4 are 
compatible with observations. Models with smaller 
lead to a decreasing initial mass with system age and are 
excluded as well. 

Our model employs a number of other parameters, such 
as the radial position and extent of the debris ring, the 
dynamical excitation and break-up strength of the parti- 
cles. We have studied the robustness of our results when 
these parameters are varied (Fig. [7|)- As long as the 
values for these parameters remain constant over age, 
varying them do not affect our conclusion on . The as- 
sumption that the dynamical excitation is constant over 
age is suspicious, in light of results from coagulation mod- 
els showing that stirring by large plantesemals increases 
gradually eccentricities of the disk particles. This is dis- 
cussed in ij5] 

There is significant uncertainty in our conclusion due 
to the small sample size. However, we argue that a larger 
sample may still not favor models with, e.g., q^ = 5. If 
93 = 5 (lower-right panel in Fig. the system that 
remains easily detectable at 2 Gyrs of age requires an 
initial solid mass of ~ 2OOM0 ~ IMj in the planetes- 
imal belt. The initial gas mass in such a belt will be 
higher than the total disk mass of a typical T-Tauri star 
(0.01 Afo).^ 

By focusing on dust luminosities, we are sensitive only 
to bodies that lie below the break-size. As seen in Fig. 

break-size marches up to few tens to a hundred kilo- 
meters by the end of a few billion years, if the disk has 
a mass of Mq = 12.4Me. 

5. DISCUSSIONS 
5.1. Coagulation Models vs. Debris Disks 

In our exercise, we have assumed a simple initial size 
distribution (eq. [T|), with all bodies larger than a few 
hundred meters described by a single power law index 
qs. We relax this assumption here. 

Simulations of planetesimal coagulation produce typi- 
cally more com plicated size distributions. For example, 
IKenvon fc Bromlev, (,2008.) started their simulations with 



all bodies at < 1 km. After tens of millions of years of 
growth, most of the mass still remains at or below 1 km, 
with only ~ 8% of the mass being accreted into bodies 
10 km or larger, ^ 6% into bodies 100 kms or larger, and 
~ 3% into bodies of order 1000 kms. We use a broken 
power-law to replicate this kind of primordial spectrum. 
We set 53 = 5.5 from 1 km to 10 k m, and q^ = 4 from 
10 km to 1000 km. Motivated by ISchlichting fc Saril 
(I20T1I) . we also consider a slightly different initial distri- 
bution with 93 = 7 from 1 km to 10 km, and 93 = 4 from 
10 km to 1000 km. Both sets of size spectrum deposit 
mass mostly at the low end (< 1 km) and little at the 
large sizes. As expected, when initial masses are deter- 
mined for different systems (Fig. IS]), we find that young 
systems require exceedingly low initial masses, while old 
systems require unphysically large initial masses. 

If we follow the luminosity evolution of such a disk, we 
will see that the disk flares brightly in the flrst tens of mil- 
lions of years, due to the large mass reservoir at the 1 km- 
range. Then the luminosity decays as (as expected 

of a (73 —4 spectrum) but with a low normalization - 
most of the disk mass has been ground down in the early 
stage and we are now left with but a scrap remnant of 
the original. Conglomeration simulations typically find 
that only a small fraction of the mass can be accreted 
to make large bo dies, before viscous stirri ng effectively 
stalls the growth. ISchlichting &: Saril ()2011|) showed that 
the fraction in large bodies can only be of order 10"'^ in 
Kuiper-belt-like environments. 

Does results in Fig. [5] allow us to exclude current con- 
glomeration models? One possible caveat in our analysis 
is the eccentricity. We discuss this below. 

5.2. Eccentricity 

We assume a static, high eccentricity (e = 0.1) for all 
systems at all times. In realistic systems, eccentricities 
can be a function of time. 

One possible cause of eccentricity evolution is coUi- 
sional cooling. 

Collisions dissipate energy, so coUisional products have 
in average lower velocity dispersion than their parent 
bodies. In a single collision, two bodies with masses mi 
and m'l (assume mi 3> m'l) impact with typical veloci- 
tie£l 

vi r^v[ ei Vkep, (8) 

where the subscript 1 indicates that this is a flrst genera- 
tion collision in our counting. Assuming that all collision 
debris fly away from the collision site with the velocity 
of the center-of-mass, i.e., all relative velocities in the 
center-of-mass frame is dissipated during the collision, 
coUisional cooling can be expressed as 



V2 



mi + m'l 



VI 1 



2mi 



(9) 



So the closer in mass the two colliding bodies are, the 
more cooling their debris experiences. If cooling domi- 
nates the eccentricity evolution, we flnd that a minimum 
eccentricity of ei = 0.13 is required (for the hard strength 
law, and ei = 0.02 for the medium law) to allow coUi- 
sional cascade to proceed all the way to micron range. 

Velocities here refer to the random component. 
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Fig. 9. — Same as the top-right panel in Fig. |6]but instead of 
a constant eccentricity (e = 0.1), here we assume that the eccen- 
tricity rises as e{t) 0.1(t/10''yrs)^''*. There is little difference 
to the inferred mass, and if anything, the data seems to argue 
that a ga = 4 model slightly overestimate the value of 53. So our 
conclusion that qs G [3.5,4] remains unchanged even considering 
eccentricity growth. 

However, even if coUisional cooling is severe, we argue 
that viscous stirring by large planetcsimals dominates 
the eccentricity evolution. This is able to raise the ec- 
centricity of collisional debris to values comparable to 
that of their parents in a time shorter than a collisional 
time. So the condition for a successful collisional cas- 
cade is reduced to e > 0.05 for the hard strength law 
and e > 0.01 for the medium strength law, i.e., the min- 
imum random motion necessary to break up the hardest 
grains (the smallest ones). 

In fact, stirring is likely to gradually raise the eccen- 
tricity of all bodies. Stirring by large bodies in the disk 
goes as e o c t^^^ (c.f. iGoldreich et al.ll2004f ). In the sim- 
ulations of lKenvon fc BromlevI (|200'8f ). planetcsimals are 
continuously stirred by Pluto-like bodies, but they only 
reach e ^ 0.1 at about a GyrsQ Under such a scenario, 
the inferred initial disk mass is similar to the original 
result (Fig. but it is clear that we prefer the same 
range of values for . 

6. SUMMARY 

Using an ensemble of bright debris disks around Sun- 
like stars, we have measured the size spectrum of their 
embedded planetcsimals. We parametrize the size spec- 
trum as dn/ds oc s""^^ and find ~ 3.5 — 4, where 
(73 = 4 corresponds to equal mass per logarithmic decade. 
The planetesimal sizes our technique probes lie between 
a couple kms to ^ 100 km. 



While this size spectrum appears consistent with re- 
sults of coagulation simulations (gs ~ 4), there are two 
lines of evidences that suggest problems in current coag- 
ulation models. 

The first line of evidence is related to the inferrc disk 
mass. The inferred initial masses for these bright disks 
are surprisinglyhigh. We find total masses reaching as 
high as lOAfgO This is comparable to the total solid 
mass in the Kuiper belt region of Minimum Mass Solar 
Nebula model, and about a factor of 100 higher than 
the mass in large Kuiper belt objects. Current coagula- 
tion models require an MMSN-like total mass to produce 
the observed density of large Kuiper belt objects. If the 
same inefficiency persists for our disks, one would require 
a total disk mass of '--^ 100 MMSN to produce those em- 
bedded planetcsimals. This is difficult to imagine. 

The second line of evidence regards the size spectrum. 
We experiment with size distributions that arise from 
coagulation simulations. We find that these distributions 
could not reproduce the luminosity distribution of the 
observed disks. Current coagulation models are highly 
inefficient in making large planetcsimals. So most of the 
mass remains at where they started, presumably ~ 1 km. 
This leads to debris disks that are too bright at early 
times and that are too dim at late times, by a couple 
orders of magnitude. 

We do not believe these discrepancies can be resolved 
by relaxing some of our model assumptions. In particu- 
lar, we argue that our estimate for (73 is unchanged even 
taking into account the fact that disk eccentricity may 
rise with time. Our results are also insensitive to the 
width of the debris ring, to the strength of bodies, and 
to the assumed upper and lower sizes. 

Because we restrict our attention to the upper envelope 
of inferred masses, our result is dominated by a handful 
of systems. Our analysis may be vulnerable to errors. 
However, the evidence is solid that debris disks remain 
fairly bright even at a few billion years. This alone dic- 
tates that there ought to be lots of mass stored in large 
(10-100 kms) planetcsimals. We address how this is ac- 
complished by revisiting coagulation model in an upcom- 
ing publication. 
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APPENDIX 



EVOLUTION OF DEBRIS DISK PROPERTIES: ANALYTICAL MODEL 

In the following, we present a simple analytical model that describes the time evolution of d ust luminosity and size 
distribution in debris disks. This model is very similar to that described in lLohne et al.l ()2008[ ). except for our choice 
for the size of the largest planetcsimals. In the following, we present results with arbitrary strength law and initial 
size distribution, followed by numerical evaluations using the hard strength law and for q^ = 4. 



^ This constraint can be reduced by a factor of unity when r a- 
diation pressure on small grains are considered HThcbaultil2009h . 

* An eccentricity of e ~ 0.3 at 40 AU corresponds to the surface 
escape velocity of Pluto. Planetcsimals have to have a near-surface 



encounter before they can reach such a high eccentricity. This takes 
time. 

^ This is for qs = 4, and even higher values are required if 
IJ3 = 3.5. 
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We approximate the body strength (eq. [S]) by two broken power-laws, 

s < Si 



{D 




(Al) 

s > Sl 



where Si is the size at which the two expressions meet. The body strength is dominated by material strength below si 



and by self-gravity above si. For the hard strength law that we adopt, a = —0.3,(3 = 1.5 and si = cm « 300 



meters. 

Combined with equation the minimum size of an impactor that causes catastrophic disruption is 

, (its) ■<■' 

^impactor — \ l \ 1+— \-^^ J 

Here, ki = 0.9 and K2 = 1-5 for our adopted strength law. 

We define a break-size, S2 = S2(t), to be the size at which the time- integrated chance of destruction per body is 
unity, or the optical depth for size S2 to be hit is, 

r{s2) = '-f. (A3) 

Bodies larger than S2 have hardly collided and they retain their primordial size distribution, while bodies smaller than 
S2 have collided many times, and they satisfy the size distribution for collisional equilibrium. If S2 > si, we adopt a 
size distribution that is piece-wise continuous. 



dn 
ds 



7i2s"«^ = n3s|""«^s-«2 si<s<S2 (A4) 

Jl^S^''^ S > S2 



where qi and 52 are the power indexes at collisional equilibrium. They are 3.5 (jDohnanvilll969n if the size ratio between 
the impactor and the target is constant. Given equation ([5]), we have qi = 3.6 and q2 = 3.0 for the hard strength law. 
This piece-wise size distribution breaks down near the blow-out size due to an abrupt deficit of small bullets. A more 
accurate derivatio n for the size distribution can be obtained by assuming that the mass loss rate is constant with size, 
as is carried out in lStrubbe fc Chiangj (|2006f ). The size distribution shows a flare-up toward the blow-out size, and the 
magnitude of the fl are- up depends on, among other things, the value of eccentricity. Our analytical results obtained 
based on equation (jA4|) should be regarded as illustrative. 

We first obtain the evolution of S2 with time. When S2 < si, i.e., collisions involve only bodies bound by the material 
strength, optical depth for destruction at S2 is determined by integrating over all its possible bullets, 

znaAa Jk^s'^^ zaAa qi — 1 

Substituting this into the definition for S2 (eq. IA3[) . we obtain 

S2 OC t -<!i+<!3-2-~i+tigi (A6) 

This yields S2 oc t^'^ for our parameters. 
Once S2 > Sl, we perform the same exercise and obtain, 

S2 OC t -<!2+<!.3-2-~2 + '»292 (A7) 

or S2 OC t"'^ for our parameters. So at early times, the break size rises steeply with time, due to an abundance of small 
bullets; while at late times, the break size rises with time more gradually due to the relative paucity of bullets. These 
two scaling relations are observed in our numerical results (Figure [5]). 

Now we proceed to derive the scaling of disk luminosity with system age. We let the infrared luminosity to be that 
portion of the starlight that is intercepted by debris particles. This is directly related to the total surface area of all 
particles, which is mostly contributed by particles around Smin0 The fractional luminosity is therefore, 

® The upper bound of the integration is chosen to be so but it . „ ■ . 

° ^ IS 01 no importance. 
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So the evolution of luminosity is dictated by the evolution of S2 with time. In particular, at late times (when S2 > si), 
the fractional luminosity decays with time gradually, 

— OC Sf~'^^ OC i5l2-<!3+2+'2-»2-!2 . (A9) 

Again, for our choice of parameters, Lir/Lh. oc If alone is varied, Lm/L^, oc t(93-3)/(2-q3) scales as 

t~-^/^ ,t~^^'^ and i~2/'^ for = 3.5,4, and 5 respectively. This forms the basis on which we decipher the primordial 

distribution of planetesimals. 

To unde rsta nd the dependence of the fractional luminosity on a range of parameters, we return to equations (jA8[) . 
T7| and (|A5[) . retaining all the neglected constants and obtaining the following expression, 



"3 

"l7 



2QAa f„ _ 1 A torb / ^_2A_ 



4a^(gi-3)''l 



„qi-<]2 



2aAa 

"3 



(-^^ - 1) ^ (ryifj; 



2 + /^i-/^igi+gi-q3 

3— ai / 
■Smin S2 < Si 

92—93 

?(92 — 1)1 2 + 92-93 + f=2-'=292 „ 

3 — gi , 

^min 52 > Si 



(AID) 



Substituting our nominal values for the indexes (ki = 0.9, K2 = 1.5, 51 = 3.6, q2 = 3.0, qg, = 4.0), we simplify the 
dependency for luminosity into (for at late times when S2 > si). 



T /A \°'^ 

„ +-0.5 fl,r0.5„-3.6 ( \ „- 1 ii,r S /( - S „-0.6 

L 



OC t-'^ '^Afo^-^a--^ '^ e-3M^B3 A-BsJ„^ (All) 



where Mq is the total mass of the disk, a its radius, Aa/a its fractional width, e the eccentricity of particles, 
the central stellar mass, s^in the blow-out size, and A, B the strengths. This relation illuminates how our procedure, 
using luminosity to infer Mq, can be affected by various parameters. For example, the actual position of the belt is a 
piece of essential information, while other values should be known roughly to within a factor of a few to avoid gross 
mis-estimate. 

Equation (jA10|) can also be used to illustrate the effect of a time- varying eccentricity on our estimate for q^. At a 
given dust luminosity, the inferred initial mass scales with the system age and the eccentricity as 

Mo oct«^-3e^('?^-3)^ (A12) 

Let the plantesimals be stirred with a time-dependence of e oc . We define a (73 as the value of q^ one obtains by 
taking a constant eccentricity (in which case Mq oc t^^~^). The true 93 is related to it as 

93 = 3 + (g-3 - 3) (^1 + ^7) • (A13) 

So for 7 = 1/4, g3 = 4, we get the true 93 = 3.75. 

Numerically we find a weaker dependency on 7. This is related to the afore-mentioned flare-up near the blow-out 
size. If instead of equation (|A4[) . we make the simplifying assumption that the mass loss rate is the same at blow-out 
size as at other sizes, but that the micron grains are destroyed by similar grains (as opposed to smaller ones), we find 
that the dust luminosity is proportional to the total number of blow-out grains, while the mass loss rate is proportional 
to the square of this number. As a result, we write 

^ oc Ml/2. (A^4) 

From this, we derive the dependence of dust luminosity on time and on eccentricity that are slightly different from those 
presented in equations (jA9p . (jAlip . (jA12[) and (jAl3p . For instance, in contrast to equation (jA13[) . the dependence of 
<73 on 7 is logarithmic. 
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